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The molecular mechanisms of coevolution between plants and insects 
remain elusive. GABA receptors are targets of many neurotoxic terpenoids, 


which represent the most diverse array of natural products known. 

Over deep evolutionary time, as plant terpene synthases diversified in 
plants, so did plant terpenoid defence repertoires. Here we show that 
herbivorous insects and their predators evolved convergent amino acid 
changing substitutions in duplicated copies of the Resistance to dieldrin 
(Rdl) genethat encodes the GABA receptor, and that the evolution of 
duplicated Rd/ and terpenoid-resistant GABA receptors is associated with 
the diversification of moths and butterflies. These same substitutions also 
evolved in pests exposed to synthetic insecticides that target the GABA 
receptor. We used in vivo genome editing in Drosophila melanogaster to 
evaluate the fitness effects of each putative resistance mutation and found 
that pleiotropy both facilitates and constrains the evolution of GABA 
receptor resistance. The same genetic changes that confer resistance to 
terpenoids across 300 Myr of insect evolution have re-evolved in response 
to synthetic analogues over one human lifespan. 


Herbivorous insects and their host plants represent the most abundant 
and diverse forms of macroscopic life on Earth. Their extraordinary con- 
temporary diversity is hypothesized to have arisen from antagonistic 
coevolution between the two wherein defence and counter-defence 
reciprocally drove diversification'?. The escape-and-radiate model’ 
of coevolution provides one historical scenario for explaining 
patterns of evolutionary diversification in plants and herbivores. It 
proposes that novel plant chemical defences evolved in response 
to herbivore-induced selection, promoting an adaptive radiation in 
a focal plant lineage. This in turn drove the evolution of new insect 
counter-defence mechanisms and an associated radiation in a special- 
ized herbivore lineage that colonizes these plants. For example, the 
ancestral Brassicales (mustards and relatives) produced glucosinolates 


asanovel chemical defence, the ‘mustard oil bomb’, through the dupli- 
cation and neofunctionalization of genes in the cyanogenic glucoside 
pathway. This was followed by the evolution of a novel detoxification 
mechanism in ancestral Pierinae butterfly larvae that permitted a host 
plant shift from Fabales to Brassicales and resulted in the adaptive radia- 
tion of the herbivores**. While consistent with the escape-and-radiate 
model, our general understanding of the genetic and molecular bases 
of plant-herbivore coevolution and their potential role in adaptive 
radiations‘ remains limited. 

Terpenoidsarethe most chemically and structurally diverse com- 
pounds knownin plants, representing approximately 60% of all known 
natural products’. The effects of terpenoids on herbivores vary from ben- 
eficial to lethal, but many terpenoids are repellent and neurotoxic’. For 
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example, naturally occurring pyrethrins are sesquiterpenoids produced 
by chrysanthemums, the powder of which has been used as a botanical 
insecticide for thousands of years’. These terpenoids impair neuronal 
function by binding to voltage-gated sodium channels. Picrotoxin, 
another natural sesquiterpenoid, acts as anon-competitive antago- 
nist (NCA) on GABA, receptors. These pentameric ligand-gated ion 
channels conduct bicarbonate and chlorideions in the central nervous 
system of vertebrates and invertebrates. In addition, a wide range of 
terpenoids act as NCAs or positive allosteric modulators (PAMs) at 
GABA, receptors, including the monoterpenoids o-thujone and thymol, 
the sesquiterpenoids bilobalide and ginkgolides, and the diterpenoids 
isopimaric acid and miltirone'??, Thus, a diversity of evolutionary 
strategies exists, through which herbivorous insects havespecialized on 
defensiveterpenoid-containing tissues from across the entire diversity 
of plants, but such counter-defence mechanisms remain poorly known. 

Unliketheir mammalian homologues in which different classes of 
GABA, subunits form a diverse family of hetero-oligomers, the struc- 
ture and assembly of insect ionotropic GABA receptors have notbeen 
fully elucidated’. The insect GABA receptor studied most intensively 
is encoded by Rdl, which was initially characterized in Drosophila 
melanogaster through a genetic screen for mutations associated with 
dieldrin resistance". Dieldrin is a synthetic organochloride (cyclo- 
diene) insecticide first produced in the USA in 1948. The Rdl-encoded 
subunits form functional homo-oligomeric chloride channels, the 
pharmacology of which is similar to that of native GABA receptors in 
insect nervous systems'*^ Other ionotropic receptor subunits, such 
as GRD (GABA and glycine receptor-like subunit from Drosophila) and 
LCCH3 (ligand-gated chloride channel homologue 3), form GABA-gated 
cation channels when heterologously expressed”. Furthermore, Ral 
orthologues have been identified in the genomes of many other insect 
species, most of which contain one copy. However, genomes of three 
moth??? and two aphid”” lineages each contain two copies of Rdl. 

Interestingly, a point mutation at position 302 that resulted in a 
single non-synonymous replacement of an alanine with serine (A302S 
or A2S, index number for the M2 membrane-spanning region) in the Rdl 
gene of some wild strains of D. melanogaster confers resistance not only 
to NCA insecticides such as dieldrin and fipronil™® but also to diverse 
plant-produced natural terpenoids, including picrotoxin, a-thujone, 
bilobalide and ginkgolides**™*. This cross-resistance pattern led us to 
hypothesize that Rdl duplications and parallel amino acid substitu- 
tions at position 2'in RDL may be a window into an uncharacterized 
resistance mechanism that evolved over deep time against neurotoxic 
terpenoids via target site insensitivity at insect GABA receptors. Here 
we addressed this hypothesis using molecular evolution, structure- 
function, diversification rate and precision genome editing approaches 
to study the effects of the resistance mutations in vivo. 


Molecular evolution of Rdlin insects 
After mining existing databases, we aligned Rdl sequences 
from 22 orders and 171 insect families together with species with 


known pesticide resistance evolution, to directly assess duplication 
and candidate site evolution. The Rdl gene is present as a single copy 
in most lineages (22 orders and 132 families). However, 2-3 copies are 
found across genomes of 52 aphid (Aphidomorpha: Phylloxeroidea 
and Aphidoidea), scale insect (Coccoidea), treehopper and leaf- 
hopper (Membracoidea), ladybird (Coccinellidae), and moth and 
butterfly (Lepidoptera) species. After the initial duplication events, 
one paralogue retaining the ancestral Rdl sequences without resist- 
ance substitutions was lost in both scale insects and ladybirds (Sup- 
plementary Table 1). In addition, in 79 transcriptomes of species 
from the above-mentioned families, we observed the same gene 
duplications and losses in 24 species; however, one ortwo paralogues 
were not found in 55 species (Supplementary Table 1). The transcrip- 
tomes from these species may be incomplete due to sequencing 
depth or sampling bias, because Rdl is primarily expressed in the 
central nervous system. Despite this uncertainty, at least eight inde- 
pendent duplications and losses of Rdl preceded the diversification 
of these insect lineages (Extended Data Figs. 1and 2, and Supplemen- 
tary Table 2). 

We found that position 2’experienced amino acid substitutions in 
all duplicated Rdl copies. 2’S was identified in Rdl2 copies of all species 
except treehoppers and leafhoppers, in which a unique substitution 
2'N was identified, and 2/Q was identified in Rd/3 copies of ladybirds 
and moths. There are double amino acid substitutions (2’P-61) in Rdl3 
copies of scale insects (Fig. 1a). 

These substitutions are exclusively associated with ancient dupli- 
cation events of Rdl. Remarkably, segregating mutations at position 
2’, such as A2'S, A2/'N and A2/G, have also been found in contemporary 
populations of more than 24 species from 6 arthropod orders in the 
past 30 yr? (Fig. 1a and Supplementary Table 3). These Rdl mutations 
confer resistance to synthetic cyclodiene and phenylpyrazole insec- 
ticides via target site insensitivity". In addition, four positions in 
MI (1276 and G279) and M3 (V339 and A343) have also experienced 
protein coding mutations in Rdl copies of different lineages, including 
treehoppers and leafhoppers, the fruitworm beetle Byturus ochraceus, 
ladybirds and moths (Fig. 1a). 

The binding site for NCAs, including picrotoxin and the syn- 
thetic cyclodiene and phenylpyrazole insecticides, is formed by the 
pore-lining 2’, 6' and 9’ residues in both insect and mammalian GABA 
receptors ™ 667 (Fig, Ib-d). Point mutations at these sites can reduce 
or abolish the potentiation effect of PAMs such as benzodiazepines, 
general anaesthetics and thymol** °°. However, the binding mode of 
benzodiazepines, such as diazepam and alprazolam, is in the extra- 
cellular domain and M3-Ml interfaces, and general anaesthetics such 
as barbiturates and propofol target the M3-MI interfaces” ”, sug- 
gesting that the substitutions in the MI and M3 regions of Rdl copies 
may be involved in PAM binding. Notably, 72 terpenoids with diverse 
structures distributed across gymnosperms and angiosperms act on 
GABA receptors as NCAs or PAMs (Supplementary Tables 4 and 5) and 
many others probably share the same modes of action. 


Fig. 1| Molecular evolution of Rdl in insects. a, Maximum-likelihood phylogeny of 
RDLusing an amino acid alignment of translated Rdl genes from resistant and non- 
resistant insect species, including Rdl duplications and amino acid substitutions 
(see Supplementary Tables 1 and 3 for all examined species). Species include those 
evolving resistance substitutions in the context of coevolutionary adaptation to 
host plants over deep time (not underlined) and those that have evolved resistance 
mutations on contemporary timescales in response to insecticides (underlined, 
n=9 species). The first set of duplicated copies are named Rdl2 (for example, 25 
and2'N) and the second set of duplicated copies are named Rdl3 (for example, 20 
and2P). Black thick branches represent inferred duplications of Rdl (medium thick 
Rdl2Zand thick Rdl3). The amino acid substitutions of six sites (1276, G279, A2’, T6’, 
V339 and A343) are highlighted in blue. Amino acid substitutions at positions 1276 
and G279 were found in the Rdl2 copies in treehopper and leafhopper species; the 
amino acid A2’was replaced in all the species containing duplicated Rdl copies 


andasubset of species that are resistant to synthetic insecticides; the amino acid 
G279 was replaced in some ladybird lineages; the amino acid V339 was replaced in 
Byturus ochraceus and ladybirds; the amino acid of A343 was replaced only in the 
Rdl3 copies of moths. Aphidomorpha, aphids (Phylloxeroidea and Aphidoidea); 
Coccoidea, scale insects; Coccinellidae, ladybirds; and Membracoidea, treehoppers 
andleafhoppers. b-d, Model ofthe picrotoxin-bound D. melanogaster RDL (insect 
GABA receptor) homo-multimer. View of the binding pocket from the parallel to 

the membrane plane (b), and side-on (c) and down-top (d) views ofthe picrotoxin- 
bound channel pore. Picrotoxinin, the main active ingredient of picrotoxin, is shown 
in ball-and-stick, side chains of amino acid residues are shown as sticks, dashed 

lines indicate hydrogen bonds, and membrane-spanning segments are indicated 

by dark-green ribbons. e f, Plot ofthe pore radii (e) and pore radius at position 2'(f) 
in wild-type and mutant RDL models. A, wild-type; S, A2’S; Q, A2/0; N, A2'N; and 

PI, A2P-T61. The dashed circle indicates the binding pocket of picrotoxin. 
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We hypothesized that the ancient amino acid substitutions at 
these positions may confer resistance to host defensive terpenoids in 
the salient herbivorous lineages. To address this, we generated models 
containing amino acid replacements docked with picrotoxin and thy- 
mol using the crystal structure of the picrotoxin and propofol-bound 
human GABA, receptor ^^! (Fig. 1b and Extended Data Fig. 3a). Our dock- 
ing simulations suggest that the 2’, 6’and 9’ M2 residues are involved 
in picrotoxin binding (Fig. 1c,d), consistent with mutagenesis, elec- 
trophysiology and structural pharmacology results'9?97?', Picrotoxin 
forms hydrogen bonds with T6'and hydrophobic interactions with 
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Locusta migratoria (Orthoptera) 
Blattella germanica (Blattodea) 
Blattella germanica (Blattodea) R 
Bemisia tabaci 

Bemisia tabaci R 

Diaphorina citri 

Daktulosphaira vitifoliae Rdl1 
Myzus persicae Rdl1 

Myzus persicae Rdl1 R 

Aphis glycines Rdl1 
Daktulosphaira vitifoliae Rdl2 
Myzus persicae Rdl2 

Aphis glycines Rdl2 
Phenacoccus solenopsis Rdl2 
Coccus sp. Rdl2 

Unaspis euonymi Rdl2 
Phenacoccus solenopsis Rdl3 
Coccus sp. Rdl3 

Unaspis euonymi Rdl3 
Nilaparvata lugens 

Nilaparvata lugens R 
Laodelphax striatellus 
Laodelphax striatellus R 
Philaenus spumarius 
Homalodisca vitripennis Rdl1 
Amrasca biguttula Rdl1 
Euscelidius variegatus Rdl1 
Homalodisca vitripennis Rdl2 
Amrasca biguttula Rdl2 
Euscelidius variegatus Rdl2 
Agrilus planipennis 

Byturus ochraceus 

Serangium japonicum Rdl2 
Harmonia axyridis Rdl2 
Coccinella septempunctata Rdl2 
Serangium japonicum Rdl3 
Harmonia axyridis Rdl3 
Coccinella septempunctata Rdl3 
Tribolium castaneum 

Tribolium castaneum R 
Acanthopteroctetes unifascia 
Neopseustis meyricki 

Plutella xylostella Rdl1 

Plutella xylostella Rdl1 R1 
Plutella xylostella Rdl1 R2 
Danaus plexippus Rdl1 

Plodia interpunctella Rdl1 
Spodoptera exigua Rdl1 
Operophtera brumata Rdl1 
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Drosophila melanogaster 
Drosophila melanogaster R 
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M1 


A2’, T6’and L9’. In M1 and M3 regions, molecular docking shows that 
the 1276 can bind thymol and other positions near the thymol-binding 
pocket (Extended Data Fig. 3b,c). The A2’S mutation was identified in 
populations of many cyclodiene-resistant insect species? and may 
result in steric hindrance of NCAs in the channel pore'^(2.3 A; Fig. 1e,f 
and Extended Data Fig. 3e) compared to the wild-type RDL (2.8 A; 
Fig. 1e, f). The2/Q and 2/Nsubstitutions further reduced the pore radius 
to 0.6 Aand0.9 A, respectively, causing more steric hindrance (Fig. 1e,f 
and Extended Data Fig. 3f,g). However, the double substitution 2’P-611 
decreased the channel diameter at position 6’ (Fig. le,f and Extended 
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Data Fig. 3h). Notably, the substitutions in M1 and M3 did not affect the 
channel-pore radius (Extended Data Fig. 3d,i-1). These results indicate 
that the substitutions could prevent the binding of the receptor with 
NCA and PAM terpenoids and may reduce the potentiation effect of 
PAM terpenoids. Overall, duplications and amino acid substitutions 
in Rdl may contribute to adaptation to insecticides and neurotoxic 
terpenoids. 


Evolution of Rdl is associated with herbivore host 
shifts and diversification 

The most ancient Rdl duplication events probably arose -298 million 
years ago (Ma) in the two common ancestors of aphids and scale 
insects, respectively (Fig. 2a). Inthe ancestral lineage of scale insects, 
two rounds of gene duplications and the loss of the ancestral version 
of the gene occurred. The Rdl gene in the ladybird lineage, which feed 
primarily on aphids and scale insects, experienced a parallel set of 
duplications and losses with respect to its prey -186 Ma. The sister 
grouptotherestoftheladybirds, the subfamily Sticholotidinae, carries 
two paralogues and has also lost the ancestral version of Rdl (Supple- 
mentary Table 1). The duplications found in treehoppers and leafhop- 
pers appeared around the same time -187 Ma or earlier since the Rdl 
sequence was not found in members of the sister family Myerslopiidae. 
Notably, the duplications that generated the Rd/2 copy are older than 
the duplications that generated the Rd/3 copy observed in the lepidop- 
teran lineages (-230 Maand -94 Ma, respectively) (Fig. 2a). Thus, these 
results suggest that most duplication events occurred near the origin 
of angiosperms”, and only duplications that generated Rdl3 in moths 
occurred after the origin of core eudicots. The relative sequential timing 
of key Rdl duplication events and subsequent terpenoid-resistance 
substitutions may dovetail with key events in the evolution of plants. 

Accordingly, we investigated in more detail the potential for coevo- 
lution between insect and plant lineages inthe context ofthe evolution 
ofterpenoid defences and Rdl evolution. In plants, terpene synthases 
are enzymes necessary for the synthesis of terpenoids. The genes 
encodingterpene synthases originated in non-seed plants -450 Ma and 
subsequently diversified in gymnosperms and angiosperms, particu- 
larly core eudicots, which together account for the extreme diversity 
of terpenoid chemicals in plants****. 

Many gymnosperm- and angiosperm-specific terpenoids serve as 
chemical defences against herbivores; for example, bilobalide and gink- 
golides in gymnosperms, and picrotoxin and anisatin in angiosperms". 
We observed that the lineages of aphids and scale insects feed far more 
extensively on gymnosperms that produce terpenoids targeting GABA 
receptors compared with related lineages. Although the hosts of most 
treehopper and leafhopper families are uncertain, species from the 
Cicadellidae also feed on gymnosperms that produce terpenoids tar- 
geting GABA receptors, while species from other related clades do not 
(Fig. 2a and Supplementary Table 6). In tandem, lepidopteran species 
with duplications of Rdl diversified on host plants that produce even more 
diverse terpenoids targeting GABA receptors (Fig. 2a and Supplementary 
Table 7). Previous studies have indicated that host plant chemistry played 
a significant role in the evolution of host shifts in herbivorous insects 
over deep time^*. Our maximum-likelihood ancestral-state reconstruc- 
tions of host preference suggest such successive host shifts from early 


non-vascular land plants (bryophytes), to gymnosperms and then to 
angiosperms over time. Each of these host plant shifts is coincident with 
the observed duplications and salient amino acid substitutions of Rdl 
(Fig. 2a). Furthermore, we found that the specific substitutions at the 2’ 
residues were strongly associated with host shifts (Extended Data Fig. 4). 

We then studied the potential macroevolutionary consequences 
of the Rdl duplications and convergent amino acid substitutions for the 
insects in which they evolved. We first applied sister-clade analysis to 
test whether ancestral shifts in Rdl copy number are associated with 
shifts in diversification rate given extant species richness. We founda 
significant increase in diversification rate using six sister-taxon pairs 
associated with transitions in Rdl copy number (Table 1 and Supple- 
mentary Table 8). 

To complement the sister-clade analysis, we then estimated net 
diversification rates across all lineages in which we were able to curate 
Rdl sequences and for subsets of lineages representing the major 
insect orders with Rd non-pesticide resistance evolution (Coleoptera, 
Hemiptera and Lepidoptera) under different extinction rate scenarios 
asin ref. 36 (096 extinction, 5096 extinction relative to speciation rate, 
9096 extinction relative to speciation rate) (Fig. 2 and Extended Data 
Fig. 2). For the full dataset, we found that the Rdl copy number 
increased the net diversification rate for each extinction scenario 
(***P < 0.001). However, this effect was primarily driven by the 
Lepidoptera (***P < 0.0001), while there was no effect inthe Coleoptera 
and a slightly negative effect on the diversification rate in the 
Hemiptera (*P « 0.05) (Fig. 2b and Supplementary Table 9). 

We further tested the relationship between Rdl copy number and 
net diversification rate using a phylogenetic generalized least squares 
methodto account for phylogenetic relatedness amongthetaxa studied. 
The Rdl copy number was positively correlated with the net diversifica- 
tion rate under all extinction scenarios (P< 0.004) and explained -1-396 
of the variation across all insect orders sampled after accounting for 
phylogenetic non-independence. This pattern was primarily driven 
by the Lepidoptera (P< 0.001) in which 10-14% of the variation in net 
diversification rates could be explained by the Rdl copy number alone. 
The effect of the Rdl copy number was not significant after accounting 
for phylogeny inthe Hemiptera- or Coleoptera-specific datasets (P> 0.1). 
These results were recapitulated under both low and high lambda trans- 
formations (Supplementary Tables 10 and 11). Overall, our results imply 
that the evolutionary genetics changes in Rdl arose as a counter-defence 
mechanism that spurred adaptive radiations in Lepidoptera, which is 
amongthe most diverse insect orders. The increased diversification rate 
in Lepidoptera may have occurred through host plant shifts that enabled 
avoidance oftheterpenoid toxins concurrently diversifyingthrough ter- 
pene synthase gene duplications in plants duringthe Mesozoic (Fig. 2). 

Many herbivores take up and store host plant toxins in defence 
againsttheir own natural enemies". However, itis not well understood 
how predators of insects have evolved the capacity to resist or tolerate 
plant defence metabolites sequestered within their prey”. Here we 
found that both scale insect prey and ladybird predators underwent 
a suite of parallel molecular evolutionary changes in sequential 
order: (1) duplication of Rdl; (2) nucleotide substitutions resulting 
in resistance-conferring amino acids in the encoded protein; and 
(3) losses of ancestral Rdl (Fig. 1a). Scale insects are primarily used by 


Fig. 2| Evolution of Rdlis associated with herbivore host shifts and 
diversification. a, Family-level ancestral-state reconstructions of feeding states 
estimated from the character states of extant species: Green, terpenoid-feeding 
herbivore on gymnosperms (in Hemiptera) and angiosperms, which produce 
diverse terpenoids targeting insect GABA receptors; red, terpenoid-feeding 
carnivore of scale insects and aphids; grey, non-terpenoid-feeding insect; and 
white, unknown. Maximum-likelihood reconstructions are shown as nodal pie 
graphs. The orange and black thick branches represent inferred duplications of 
Ral (medium thick Rdl2 and thick Rdl3). Family names in black indicate that the 
Rdlsequence from at least one species in the family lineage is known. Green, pink, 


and purple background colours mark Hemiptera, Coleoptera and Lepidoptera, 
respectively. The heat map shows the percentage of terpenoid-feeding species 
(Supplementary Tables 6 and 7) and net diversification rate. The timeline relates 
the convergent evolution of Rdl duplications in insect taxa to the origin of 
gymnosperms and angiosperms, and the subsequent diversification of each. The 
evolution of terpene synthase genes (TPS) in plants is taken from ref. 34. The tree 
was rooted with sequences from Dermaptera (Diplatys sp.), Plecoptera (Moselia 
infuscata), and Orthoptera (Myrmecophilus sp.) (Supplementary Table 1). 

b, Relationship between Rdl copy number and net diversification rate across all 
insects surveyed for Rdl copy number variation and resistance evolution. 
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Table 1| Sister-clade analysis results 


Species pair removed from sister-clade analysis Chi-squared d.f. P value 

None 90.03757 1 2.3368x10 7 
Aphidamorpha/Coccoidea: Psylloidea 43.84667 1 3.551349x10" 
Coccoidea: Aphidamorpha 35.41264 1 2.667503x10 ?? 
Membracoidea: Cicadoidea/Cercopoidea -8.580884 1 1 
Coccinellidae: Cerylonidae/Bothrideridae -8.581524 1 1 

Heteroneura: Lophocoronoids/Hepialoids -10.81037 1 1 
Noctuoids/Lasiocampoids/Bombycoids: Drepanoids -10.74069 1 1 


We performed a full sister-clade analysis using all species pairs (no species pair was removed from sister-clade analysis). We also performed additional sister-clade analyses by removing 
one pair of sister clades to test for the contribution of any one pair on the full analysis (see Methods for more information). P values represent a one-sided test and were calculated using the 


‘richness.yule.test’ function in the R package ape. 
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Fig. 3| Amino acid replacements and duplications of Rdl contributing to 
picrotoxin and thymol insensitivity in gene-edited D. melanogaster. 
a,b, Adult survival of flies reared on diets containing picrotoxin (a) or thymol (b) 
among distinct lines. Rdl genotype: +/+, w™° wild-type; +/+*, A2’A engineered 
control; S, A2/5; Q, A2’Q;N, A2N;P, A2P; and I, T61. log(inhibitor) versus response 


nonlinear fit was performed for a. One-way ANOVA with post hoc Bonferroni 
correction was used for b (Supplementary Table 16); mean + s.e.m., n = 12, 3,3, 
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3,5,3,6, 6,12, 6,10, 5,3,3,3,3,3, 4 and 6 biological replicates, and 10 females per 
replicate; *P < 0.05, **P< 0.01, ***P < 0.001, ****P < 0.0001. LC;o, median lethal 
concentration. c, Model of stepwise duplications of Rdl in the insect lineages 
generates multiple subunits to form candidate homo- and hetero-pentameric ion 
channels conferring resistance to neurotoxic terpenoids. The line break between 
A, Sand Qin Lepidoptera indicates that these are not on the same chromosome. 


ladybirds as prey’. Our results suggest an evolutionary scenario in 
which the duplications and amino-acid-changing substitutions of Rdl 
in the ladybird ancestor conferred resistance to terpenoids used by 
scale insects as a chemical defence against natural enemies, which in 
turn contributed to the evolutionary radiation of ladybirds (Fig. 2a). 


Amino acid substitutions of Rdl contribute to 
terpenoid insensitivity 

To investigate whether the Rdl substitutions we identified provide 
resistancein vivoto both NCA and PAM terpenoids in whole organisms, 


and to assess their potential pleiotropic effects, we introduced 
them into the native Rdl gene of D. melanogaster mainly using a two- 
step CRISPR-Cas9 genome editing approach (Extended Data Fig. 5). To 
recapitulate the effects of gene duplication events and assess pleio- 
tropic effects of different combinations of mutations, we also made 
a series of heterozygous lines to study the particular Rdl genotypes 
we found across a diversity of insect lineages (Fig. 3). 

We found that three lines (N, 61and P-61) were homozygous lethal, 
suggesting that these mutations exact high fitness costs. We then 
measured survival and performance using picrotoxin and thymol as 
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representative NCA and PAM terpenoids, respectively, in different 
feeding assays. We found that S/+ (Ser/Ala) and S/S increased adult 
survival compared with the wild type when exposed to increasing 
concentrations of picrotoxin in the diet. Q/+ (Gln/Ala) also increased 
adult survival when fed on picrotoxin but were more sensitive to 
picrotoxin than S/+ and S/S. As expected, S/Q and Q/Q increased 
adult survival compared with the others at 10 mM picrotoxin (Fig. 3a 
and Extended Data Fig. 6a). Consistent with findings from a study 
that created ‘monarch flies’ that were fully resistant to dietary card- 
enolides", the increased adult survival we observed was not due to a 
reduction in feeding rate or toxin ingestion (Extended Data Fig. 7). The 
N/+ (Asn/Ala) and S/PI (Ser/Pro-6'lle/Thr) genotypes also increased 
adult survival when fed picrotoxin (Fig. 3a and Extended Data Fig. 6b). 
In addition, we tested mutants in which Rdl genotypes (S/P, SI/+, PI/+, 
P/+,1/+ and P/P) are absent in the scale insect lineages. The homozygous 
P/P genotype increased adult survival upon picrotoxin exposure and 
the heterozygous S/P and SI/+ (S/+-611/+) both showed higher resist- 
ance than S/PI flies. However, three additional heterozygous lines, PI/+, 
P/+ and I/+, showed no significant resistance to picrotoxin (Extended 
Data Fig. 6b). In contrast, most of these lines were resistant to varying 
degrees upon exposure to thymol, except P/P, which was more sensi- 
tive than the wild type (Fig. 3b and Extended Data Fig. 6c,d). We also 
tested whether the amino acid substitutions in M1 and M3 regions 
provide resistance to thymol, but we found no evidence for an effect 
(Fig. 3b and Extended Data Fig. 6e). Our results indicate that amino acid 
replacements in the M2 of RDL encoded by the Rdl paralogues confer 
resistance to these exemplar NCA and PAM terpenoids in vivo, and 
the loss of the sensitive copy of Rdl (ancestral character state) may 
further enhance resistance. 

Coexpression of Rdl1 and Rdl2 from moths can generate 
hetero-pentameric GABA receptors with intermediate sensitivity to 
dieldrin compared with their homo-oligomers?". More importantly, a 
duplication of the Rdl locus (2A in one copy and 2/S in the other copy) 
resulted in functional and permanent heterozygosity with intermedi- 
ate insecticide resistance in two natural D. melanogaster strains?. 
In agreement with this, we found that the diverse RDL subunits may 
formhetero-oligomers that result in variation interpenoid resistance 
(Fig. 3a,b and Extended Data Fig. 6). Therefore, the duplications of 
Rdl in these insect lineages could generate a multiplicity of RDL sub- 
units that form homo- or hetero-pentameric receptors. For example, 
two Rdl copies could result in 2° (32) unique combinations and three 
copies could result in 3° (243) unique combinations, resulting in poten- 
tially diversefunctional properties and pharmacological characteris- 
tics (Fig. 3c). Our findings suggest that subtle genotypic differences in 
Rdl belie a far more complex range of phenotypic possibilities at the 
level of the pentameric receptor in this system. However, it should be 
noted that the subunit stoichiometry of any insect GABA receptor is 
unknown and the actual number of receptor subtypes could be fewer 
thantheabove calculations. Further physiological and structure-func- 
tionstudies ofthe GABA receptor are needed to evaluate the biological 
significance of these findings. 


Amino acid replacements of Rdl are associated 
with fitness costs 

Given the essential role of GABA receptors in the nervous system and the 
high degree of RDL sequence conservation across animals, mutations 
that confer resistance may impair the ancestral function. We addressed 
this by measuring potential fitness costs in the Drosophila mutants 
we created. In comparison with control flies, Q/Q and P/P genotypes 
decreased egg laying of females, while other lines showed no measured 
defects (Extended Data Fig. 8). Previous studies found that the A2'S 
mutation causes a temperature-sensitive phenotype in Drosophila ^^^, 
so we also performed a heat-shock test across our lines. We found that 
S/Sgenotypeflies beganto exhibit paralysis behaviour after exposure 
to 38°C for 15 min. Notably, Q/Q homozygote flies became paralysed 


more rapidly and were completely unable to move after 5 min (Fig. 4a 
and Supplementary Video 1). By contrast, S/+ and Q/+ heterozygous 
flies showed higher thermal resistance compared with S/S and Q/Q 
genotype flies, respectively. The S/Q heterozygote flies displayed 
an intermediate heat-shock sensitivity that was in-between the 
levels found in the S/S and Q/Q homozygous flies (Fig. 4a). The N/+ 
heterozygous flies also exhibited defects compared to wild-type 
flies (Fig. 4a). Importantly, S/PI heterozygous flies displayed higher 
temperature tolerance than S/P, P/P or P/+ flies, suggesting that the 
pleiotropic costs of 2P are rescued to some degree by the 6/1 mutation 
(Fig. 4b). Finally, we measured locomotory behaviour using an 
automated monitoring system and found that all lines carrying 
engineered mutations exhibited a decrease in locomotion, particularly 
the S/PI and P/P lines (Fig. 4c and Extended Data Fig. 9). The severe 
deleterious effects of S/PI on movement may be tolerable in the scale 
insect radiation since the immature life stages most susceptible 
to predation and parasitism are slow moving or immobile. Overall, 
these results show that the coexistence of multiple Rdl copies amelio- 
rates antagonistic pleiotropy while maintaining resistance to diverse 
terpenoids. 


Discussion 

Identifying the molecular mechanisms arising from coevolution 
between plants and herbivorous insects is difficult because it involves 
determining how sets of species have coadapted through vast expanses 
of time'^*^*^, The repeated evolution of resistance to synthetic NCA 
insecticides provides a useful window into how insects responded to 
repeated bursts of novel terpenoid insecticide production across deep 
timeasancient plants diversified and their terpene synthase repertoire 
expanded via genome duplication events^^. 

Amino acid substitutions and gene duplication events repeat- 
edly evolved in the insect GABA receptor subunit gene Rd/ through 
selection for cyclodiene and phenylpyrazole resistance over the past 
70 yr/*5? The A2'S and A2'N substitutions, which confer resistance 
to these insecticides, are now globally widespread in diverse arthro- 
pod lineages”. Remarkably, the same mutations and parallel gene 
duplication events became fixed in four different herbivorous lineages 
and one predator lineage over the past 300 Myr, indicating that the 
adaptive path is repeatable and therefore predictable under similar 
selective pressures*' 5, 

These molecular evolutionary changes in Rdl were associated 
with host plant switching events across herbivorous insects and 
an increased net diversification rate in the Lepidoptera, one of the 
most diverse insect orders (Fig. 2). However, it remains unclear why 
Coleoptera and Hemiptera lack a signal of increased diversification 
rate associated with the evolution of resistant GABA receptors. The 
Lepidoptera experienced the most recent Rdl duplications, which 
occurred -94 Ma, while Rdl duplications in the Coleoptera and 
Hemiptera occurred over 180 Ma. This, together with the long 
internal branches and likely high extinction rates in the Coleoptera 
and Hemiptera over deep time, may have obscured signals of a 
diversification rate increase associated with Rdl duplication. Alterna- 
tively, adaptive evolution of Rdl did not spur diversification in these 
two orders. 

Fossil evidence suggests that the stem lineage of aphids, which 
originated in the Triassic, fed on gymnosperm hosts^^!, Our results 
provide a conceptual model of how this may have proceeded at the 
geneticand protein receptor levels in several major herbivorous insect 
lineages as well as the ladybirds, a major predator of sternorrhynchan 
insects (aphids, scale insects, psyllids and whiteflies). Amino acid 
replacements at position 2' occurred in all duplicated copies and these 
RDLsubunits may form heteromericion channels that confer resistance 
toterpenoids while minimizing fitness costs through the amelioration 
ofantagonistic pleiotropy. We hypothesizethat such properties would 
provide persistent fitness advantages in ancient herbivore lineages that 
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Fig. 4| Amino acid replacements and duplications of Rdl are associated 

with fitness costs in D. melanogaster. a,b, The percentage of flies observed 

to be walking when exposed to high temperature for 30 min. Repeated- 
measures ANOVA (general linear model) and post hoc Bonferroni correction 
(Supplementary Table 16) for a, n = 5,5,5,5, 4, 5, 3and 5 trials; for b, n = 5,5,3,5,5, 
5,5, 5andS5trials, and 10 females per trial, mean + s.e.m. c, Daily crossing activity. 
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One-way ANOVA and post hoc Bonferroni correction (Supplementary Table 16) 
were used for statistical analysis; mean + s.e.m., n = 64, 58, 63, 62, 61, 64, 54, 

64, 62, 64, 64, 64, 64 and 59. Genotypes as in Fig. 3. Letters indicate significant 
differences among genotypes and experimental groups denoted by ‘ab’ are not 
significantly different from either a or b groups. 


encountered an increasing diversity of terpenoid toxins targeting the 
GABA receptor through both NCA and PAM mechanisms. 

Unlike other insect lineages with two Rdl copies, the macro- 
evolutionary patterns underlying terpenoid resistance in moths 
carrying multiple copies may reflect two peaks on an adaptive land- 
scape. Given that the Rdl3 copies carrying 2’Q confer higher resistance 
to picrotoxin than the 2’S-harbouring Rd/2 (Fig. 3a and Extended Data 
Fig. 6a), a stepwise evolutionary path was taken by the moth radia- 
tion. The initial duplication that generated Rd/2 may have provided an 
advantage that spurred early host shifts to gymnosperms. The origin 
of the angiosperms and the rapid diversification of core eudicots 
provided more diverse host plants, terpenoid synthases and terpe- 
noids by the early Cretaceous. However, many lepidopteran clades 
are restricted to plants from particular angiosperm lineages”. The 
second duplication event that generated Rdl3 may have facilitated 
an escape from ancestral host plant lineage of moths and exploi- 
tation of other core eudicots, which then led to a second major 
adaptive radiation. This hypothesis is congruent with patterns of 
host use and diversity of the superfamily Noctuoidea, the most 
species rich of any superfamily, containing nearly one-third of all 
lepidopteran species”. 

Core eudicots produce a greater diversity of PAMs compared with 
gymnosperms. Structures of GABA, receptors in complex with PAMs, 
suchas benzodiazepines, suggest that the binding pocket includes the 
M3-M1 interfaces along with the 15’ residue in M2 (ref. 27). However, 
we found that the substitutions in M1 and M3 have no effect on PAM 


resistance (Fig. 3b and Extended Data Fig. 6e), but the substitutions 
at the 2’ position in M2 confer resistance to the PAM thymol (Fig. 3b). 
These M2 substitutions both directly affect the binding at the NCA drug 
binding site and also allosterically destabilize the PAM drug preferred 
desensitized state as previously described”. 

Gene duplications and amino acid substitutions in Rdl, together 
with the pre-mRNA A-to-I editing and alternative splicing’, may also 
confer structural diversity to the pentameric GABA receptor com- 
plex, yet the mechanistic details of how homo- or hetero-pentameric 
receptors assemble are unclear. Such cryptic GABA receptor diversity 
would facilitate the maintenance of terpenoid resistance but also relax 
constraints imposed by resistance substitutions over deep time. 

In summary, we found evidence that the convergent evolution of 
Rdl conferred resistance to neurotoxic plant terpenoids associated 
with herbivore host shifts, higher net diversification in moths and 
butterflies, and an adaptive radiation of a predatory beetle family. 
Our results support a cascading model of coevolution in which the 
increasing complexity of terpenoid defensive cocktails in plants was 
followedby the evolution of resistant GABA receptors in the herbivores 
and their predators. Thus, terpenoids and resistant GABA receptors 
could be one example of coupled key innovations that have contributed 
to the adaptive radiation of insects and plants at macroevolutionary 
timescales. A scenario similar to the one that unfolded over 300 Myr 
of coevolution has also played out across a single human generation, 
as insects adapted to synthetic insecticides that functionally mimic 
natural neurotoxic terpenoids. 
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Methods 

Identification of Rdl genes in insects and phylogenetic 
analyses of RDL 

To identify Rdl genes in insects, we performed a two-step analysis: 
(1) we used D. melanogaster, Bombyx mori, Chilo suppressalis and 
Acyrthosiphon pisum Rdl genes as queries to perform BLASTp 
and TBLASTn searches against genomes and transcriptomes from 
GenBank, AphidBase, InsectBase 2.0, Fireflybase, DRYAD, Lepbase 
and GigaDB, because the Rdl genes from these species have already 
been confirmed???" (2) we verified the candidate genes by BLASTp 
again without a species limit as we have shown before^*^. We took 
all the candidate Rdl genes that were reciprocal best hits with the 
D. melanogaster and B. mori Rdl genes and then renamed them 
(Supplementary Table 1). The Rdl copies were mapped on chromo- 
somes in three representative species (Extended Data Fig. 10). 

To estimate the evolutionary relationships of Rdl paralogues 
among different orders, we used IQ-TREE (v.1.6.12) to build amaximum 
likelihood (ML) phylogenetic tree using an amino acid alignment of 
RDL”. To create the alignment, nucleotide sequences were predicted 
andtranslated to proteins by NCBI ORFfinder (https://www.ncbi.nlm. 
nih.gov/orffinder/). The four Varroa destructor (Arachnida) RDL amino 
acid sequences were used as the outgroup. We then selected the RDL 
sequences of more than 360 amino acids in Hemiptera, Coleoptera, 
Hymenoptera and Lepidoptera and estimated the phylogeny. We report 
the Rdl sequences used in our analysis in Supplementary Table1. Protein 
sequence alignment was then performed using MUSCLE in MEGA X”. 
The best-fit model was VT+1+G4 estimated using ModelFinder®, and an 
ML tree inference was subsequently conducted in IQ-TREE with default 
parameters. The tree with the best log-likelihood was selected as the 
tree for further analysis. 

We found that duplications of Rdlare always associated with amino 
acid substitutions at site 2’. NCA insecticides, such as dieldrin and 
fipronil, also act on this site and many agricultural pests have already 
evolved segregating mutations at position 2’, which confer resistance to 
insecticides". In light of this initial finding, we decided not to analyse the 
2'substitutions intheabsence of duplication events. In particular, we did 
not find parallel amino-acid-changing mutations in sister groups with- 
out Rdl duplications, implying that these may have arisen recently and 
are associated with insecticide resistance. For example, in Thripidae, 
we observed that Rdl sequences from Frankliniella occidentalis, 
Thrips tabaci, and T. palmi have the A2’S mutations, whereas those 
from F. cephalica and Aptinothrips rufus do not. Notably, F. occidentalis 
from Utah (USA, accession number: GCYRO1013327.1) does not have 
A2’S Rdl mutations, whereas F. occidentalis from Yangzhou has an A2’S 
mutation associated with insecticide resistance (China, accession num- 
ber: MH249048.1), suggesting that the F. occidentalis from Yangzhou 
recently evolved resistance to insecticides via target site insensitivity. 


Species tree 
To further study the duplication histories of Rdl across insects and con- 
duct phylogenetically informed diversification studies, we estimated a 
species tree of all lineages screened for Rdl copy number variation and 
resistance substitutions using 1,123 single-copy orthologues extracted 
using the BUSCO pipeline and the insecta_odb10.2019-11-20 marker 
set^?. We used genome and transcriptome assemblies downloaded 
from GenBank, AphidBase, InsectBase 2.0, Fireflybase, DRYAD, Lepbase 
and GigaDB. Amino acid sequences of each orthologue were aligned 
across all genome and transcriptome assemblies in which the ortho- 
logue was complete. We aligned amino acid sequences in MAFFT? and 
then inferred ML trees for each orthologue in IQ-TREE after automatic 
best-fit model selection?*. These gene trees were then used to infer a 
species tree in ASTRAL”. 

Afamily-level phylogeny wasthen estimated on the basis of previ- 
ously published sources (see Supplementary Table 2 for references), in 
which we pruned the taxa to a single species representing each family. 


Threeorders, Hemiptera, Coleoptera, and Lepidoptera, were used for 
thisanalysisandthetree was rooted with sequences from Dermaptera 
(Diplatys sp.), Plecoptera (Moselia infuscata), and Orthoptera (Myrme- 
cophilus sp.) (Supplementary Table 1). 


Diversification 

We applied sister-clade analysis using the ‘richness.yule.test’ func- 
tion in the R package ape to test whether ancestral shifts in Rdl copy 
number are associated with shifts in diversification rate given extant 
species richness. We used stem ages taken from primary phylogenetic 
studies (see Supplementary Table 2 for references) and species rich- 
ness estimates from the Global Biodiversity Information Facility”. 
Rdl copy number is a good proxy for insecticide resistance because 
duplications of Rdl are tightly associated with resistance substitu- 
tions. Sister-clade analysis is a standard approach to understanding 
asymmetry on a phylogenetic tree without requiring fully sampled 
phylogenies, which is useful in this case because our sampling of the 
insect phylogeny is incomplete. Although sister-clade analysis has a 
long history, biases exist in comparisons between sister clades with 
different crownages”. The reason for this bias is that younger clades are 
more likely to be species poor and contain the derived trait because the 
waiting time given by the stem branchis long. In our analysis, we used 
stem ages and model species richness probabilities as a Yule process”. 
We conducted the sister-clade analysis using six sister-clade pairs and 
tested for the contribution of any one pair by conducting additional 
sister-clade analyses with one pair removed. 

Next, we estimated net diversification rates using the ‘bd.ms’ 
function in the R package geiger*' under different extinction rate 
scenarios (0% extinction, 50% extinction relative to speciation rate, 90% 
extinction relative to speciation rate). Stem ages and species richness 
estimates were taken as above. 


Ancestral-state reconstruction of feeding states 

Terpenoids targeting GABA receptors were collated from previous stud- 
ies on the basis of pharmacological assays (Supplementary Table 4). 
The distribution of these terpenoids was primarily taken from ref. 68 
(Supplementary Table 5). For the family-level feeding preference of 
Hemiptera, we calculated the proportion and the number of species 
that feed on gymnosperms (Supplementary Table 6). Data were taken 
from ref. 69. For the family-level feeding preference of Lepidoptera, 
we calculated the proportion and the number of species that feed 
on terpenoid-containing plants (Supplementary Table 7). Data were 
taken from the ref. 70. Any family in which <3 species were excluded 
for family-level ancestral-state reconstructions and in which there is 
a >1% result indicates that there is strong evidence of feeding states. 
Thedietary preferences ofthe ladybird species were based on previous 
studies (Supplementary Table 12). 

We conducted ancestral-state reconstructions on the family-level 
phylogeny (the species tree) using PastML with default parameters". 
In addition, ML and marginal posterior probabilities approximation 
methods were used to assess the state posteriors on a node. We visua- 
lized the ancestral states and phylogenetic trees using the Interactive 
Tree Of Life”. 


Associations of RDL M2 sequences and feeding states 

We used TraitRateProp to detect the associations between genotypes 
and phenotypes across the phylogeny”. We focused on the association 
between the M2 sequences and feeding states. A rooted ultrametric tree, 
RDL amino acid alignments and trait states were used for the analysis. 
The rooted ultrametric tree was constructed on the basis of previously 
published sources (see Supplementary Table 2 for references). 


Fly strains 
Drosophila melanogaster flies were reared on conventional 
cornmeal-agar-molasses medium at 60% + 10% humidity and 25 +1°C 
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with the day and night periods of the 12 h:12 hlight:dark cycle. The light 
was set at 07:00. The w™” (#5905) strain was used as the wild-type strain 
inthis study. The following strains were obtained from the Bloomington 
Stock Center (Indiana University): Ral/** (RAIS, #1675), vas-Cas9 
(#51323, #51324), and nanos-Cas9 (ref. 74). 


Generation of knock-in flies 

To generate the Rdl knock-in lines, we used a two-step CRISPR (clustered 
regularly interspaced short palindromic repeats)-Cas9 method relying 
on homology-directed repair (HDR). We designed single guide (sg) 
RNAs using E-CRISP (http://www.e-crisp.org/E-CRISP/) and CHOPCHOP 
(https://chopchop.cbu.uib.no/)””. The sgRNAs were synthesized by 
GenScript and then subcloned into the PCFDS vector for expression. 

We used two specific approaches to generate Rdl knock-in lines 
(Extended Data Fig. 5 and Supplementary Tables 13-15). First, since 
many terpenoids, including picrotoxin and thujone, and synthetic 
insecticides, act on the 2’, 6’ and 9’ residues of M2, we reasoned that 2’ 
mutations could confer resistance to these chemicals. Therefore, we 
synthesized a single 1,000 bp homology template from D. melanogaster 
genomic DNA containing A2’Q or A2/P mutation and subcloned it with 
EcoRI and Hindlllsites into the pUC57 vector. Following embryo injec- 
tion procedures described below, the knock-in flies were identified 
through a screen in which we reared adult flies on a picrotoxin diet 
and genotyped Rdl using Sanger sequencing of PCR amplicons from 
the surviving flies. We used this method to generate the Rd/?? strain. 

Second, we synthesized two 1kb genomic DNA fragments 
obtained from D. melanogaster genomic DNA corresponding to the 
M2 sequences ofthe Rdl gene as homology arms and subcloned them 
into the pBSK-attP-3xP3-RFP-loxP vector” to generate the Rd?” 
strain. Next, the attP-3xP3-RFP-loxP sequence was replaced with each of 
the point mutation alleles in Rd/ through homologous recombination, 
generally following ref. 41. Candidate 3xP3-RFP-positive or -negative 
flies were verified by PCR of the targeted region followed by Sanger 
sequencing. We used this method to generate the Rd? ?*F*, RAMYA, 
Rd", Rdl”, Rd Tél Rdl", R dl, R dics, Rdl ZES, Rdl”, and 
Rdl?*" strains. Note that Rdl“ ^ is an engineered control strain contain- 
ing synonymous codon mutations at PAM sites, which was created to 
account for potential pleiotropic or off-target effects of the two rounds 
of CRISPR-Cas9 mediated HDR. 

To generate these knock-in lines, a plasmid mixture with the donor 
vector and sgRNAs was injected into vas-Cas9 (51323), nos-cas9 (78781) 
or RAPP2REP/yas-Cas9 (51324) embryos, as shown in Extended Data 
Fig. 5and Supplementary Table 13. The vas-Cas9 (51324) strain carry- 
ing the 3xP3-GFP marker on chromosome3 was crossed with Rd [^ FP" 
to produce heterozygous flies as Go, for embryo injections. Typically, 
250-300 embryos were injected (UniHuaii) and RFP-positive flies 
were screened under a fluorescence microscope and crossed with 
double-balanced flies (TM3/TM6B). The G,, were then screened for 
RFP- and GFP-negative flies, which were crossed with double-balanced 
flies. The G,, were verified by PCR and Sanger sequencing as above. 


Bioassays 

Picrotoxin (Tokyo Chemical Industry, C0375) was resuspended and 
serially diluted in 10096 acetone (Sinopharm, 10000418) solution. 
Thymol (Aladdin, T104426) was resuspended in 100% acetone and 
serially diluted in 50% v/v acetone/50% deionized water solution. All 
compounds were freshly prepared before each assay. Then, 100 pl of 
picrotoxin or thymol solution was added to the surface of 12 x 78 mm 
Drosophila vials containing 900 pl 2% (w/v) agarose and 5% (w/v) 
sucrose. We prepared the vials 20-24 h before the assays to allow for 
the complete distribution of the chemicals and the evaporation of 
the solution as previously described”. Ten mated female flies aged 
5-7 days old were gently introduced into the vials and reared at 25 °C, 
30% humidity and a 24 h (12 h:12 h dark:light) period. Adult survival 
was monitored for 48 h. 


Behavioural assays 

For egg-laying assays, 10 females and 5 males were grouped and reared 
on conventional cornmeal-agar-molasses medium in single vials. After 
5-7 days, the female flies were transferred to 35 mm tissue culture 
dishes containing fresh food for 24 h, and the number of eggs was 
manually counted under a stereo microscope. 

For temperature sensitivity assays, 10 female flies aged 5-7 days 
old weregently introduced into new empty glass fly vials by aspiration. 
The flies were exposed to a 38 ?C temperature and recorded with a 
video camera for 30 min. The percentage of flies awake in each tube was 
measured every 30 s. The knock-in lines exhibited paralysis behaviour 
when exposed to the high temperature, in which the animal lies on its 
dorsum and is unable to locomote. 

For locomotion assays, locomotor behaviour was performed 
as previously described™. Individual 5- to 7-day-old virgin female 
flies loaded into tubes with 2% (w/v) agarose and 5% (w/v) sucrose 
were monitored using the Drosophila Activity Monitoring System 
(DAMS, Trikinetics) in 24 h light/dark cycles. The flies were allowed 
to adapt to the new environment for 1 day and data were collected 
in 30 min bins. Average daily activities were measured on the basis 
of a2 day window. 


Molecular docking 

The human a1(3y2 and a102y2 GABA, receptors and picrotoxin and 
propofol-bound crystal structures (PDB: 6HUG and 6X3T) provided 
homology templates for building a homo-pentameric model of the 
Drosophila RDL receptor??? Molecular Operating Environments 
(MOE, 2015.10) was used as previously described**. The best model 
was chosen on the basis ofthe scoring values and evaluated using the 
UCLA-DOE server and Ramachandran plots. The models of RDL muta- 
tions were generated using Swiss-PdbViewer to introduce the amino 
acid substitutions and minimize the energy ofthe resulting structures”. 
The pore diameters of the models were calculated using the HOLE in 
WinCoot??. The MOE-Dock programme was used for docking with 
default parameters. 3D protonation was added, water molecules were 
deleted, and the energy ofthe models and ligands was minimized. The 
ligand was kept flexible during molecular docking calculations. A lower 
binding-free energy structure, which indicates a better interaction 
between the model and ligand, was chosen for the analyses. 


Statistical analysis 

All statistical analyses were performed using GraphPad Prism 7 
(GraphPad software). The adult survival curves were analysed using 
log(inhibitor) versus response nonlinear fit and thelog-rank (Mantel- 
Cox) test, and P values were compared to the Bonferroni correction 
values: Qadjustea = 0.05/(n(n-1)/2), where n is the number of groups. 
Repeated-measures analysis of variance (ANOVA) and post hoc 
Bonferroni correction were performed for temperature sensitivity 
assays. One-way ANOVA and post hoc Bonferroni correction were 
used for locomotion assays. Finally, Kruskal-Wallis and post hoc 
Mann-Whitney Utests were performed for egg-laying assays. Details 
on theother statistical methods are reported in the figure legends. 


Reporting summary 
Further information on research design is available in the Nature Port- 
folio Reporting Summary linked to this article. 


Data availability 
Net diversification rates and sequence alignments are available as 
Supplementary Data. Source data are provided with this paper. 
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Varroa destructor Rdl4 
Varroa destructor Ral2 
Ptyssoptera sp. Rall 
Plutella xylostella Rdl1 
Papilio machaon Rdl1 
Heliconius erato demophoon Rdl1 
Heliconius erato lativitta Rdl1 
Heliconius cydno Rdl1 
Danaus plexippus Rdl1 
Papilio glaucus Rdl1 
Papilio polytes Rdl1 
Phoebis sennae Rdl1 
Pieris rapae Rdl1 
Calycopis cecrops Rdl1 
Lerema accius Rdl1 

Chilo suppressalis Rdl1 
Drepana arcuata Rdl1 
Spodoptera exigua Rdl1 
Helicoverpa armigera Rdl1 
Amyelois transitella Rdl1 
Ostrinia furnacalis Rdl1 
Plodia interpunctella Rdl1 
Antheraea yamamai Rdl1 
Junonia coenia Rdl1 


p-— Bombyx mandarina Rdl1 


Bombyx mori Rdl1 
Manduca sexta Rdl1 
Helicoverpa armigera Rdl3 
Spodoptera exigua Ral3 
Bombyx mandarina Rdl3 
Bombyx mori Rdl3 
Antheraea yamamai Rdl3 
Manduca sexta Rdl3 
Ptyssoptera sp. Rdl2 
Ostrinia furnacalis Rdl2 
Papilio glaucus Rdl2 

Plodia interpunctella Rdl2 
Amyelois transitella Rdl2 
Bombyx mandarina Rdl2 
Bombyx mori Rdl2 
Spodoptera exigua Rdl2 
Helicoverpa armigera Rdl2 
Antheraea yamamai Rdl2 
Manduca sexta Rdl2 

Chilo suppressalis Rdl2 
Papilio machaon Rdl2 
Papilio xuthus Rdl2 

Papilio polytes Rdl2 

Lerema accius Rdl2 

Danaus plexippus Rdl2 
Heliconius erato demophoon Rdl2 
Heliconius erato lativitta Rdl2 
Heliconius cydno Rdl2 
Drepana arcuata Rdl2 
Junonia coenia Rdl2 
Phoebis sennae Rdl2 

Pieris rapae Rdl2 

Plutella xylostella Rdl2 
Calycopis cecrops Rdl2 
Harmonia axyridis Rdl2 
Stethorus sp. Rdl2 
Harmonia axyridis Rdl3 
Stethorus sp. Rdl3 

Biphyllus sp. Ral 

Tribolium castaneum Rdl 
Aethina tumida Rd! 
Leptinotarsa decemlineata Rdl 
Anoplophora glabripennis Rd! 
Diabrotica virgifera Ral 
Agrilus planipennis Ral 
Photinus pyralis Ral 
Nicrophorus vespilloides Rdl 
Apis mellifera Rdl 


L— Bombus impatiens Rdl 


Camponotus floridanus Rdl 
Harpegnathos saltator Rdl 
Polistes dominula Rdl 
Athalia rosae Ral 

Brasema neomexicana Ral 
Bootanomyia dorsalis Ral 
Brachymeria minuta Rd! 
Perilampus aeneus Rdl 
Orasema simulatrix Rdl 
Torymus bedeguaris Rd] 
Nasonia giraulti Ral 

Nasonia vitripennis Rdl 
Pteromalus puparum Ral 
Pachycrepoideus vindemmiae Rdl 
Trichogramma chilonis Rdl 
Trichogramma evanescens Rdl 
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Diglyphus isaea Rdl 
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Encarsia formosa Rdl 
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Extended DataFig. 1| Phylogenetic relationships of RDL in insects. Maximum subunit sequences ofthe mite Varroa destructor were used as outgroups. The 


likelihood based phylogenetic tree of RDL from four orders: Hemiptera, valueson the branches represent bootstrap support. The deep divergence of 
Hymenoptera, Coleoptera, and Lepidoptera. Protein sequence alignment Rdl duplications in three orders suggests duplications have occurred before the 
was performed using MUSCLE in MEGA X. The best-fit model was VT + I + G4 diversification of their respective lineages. For a complete list of taxa surveyed 
estimated using ModelFinder, and a maximum likelihood (ML) tree inference see Supplementary Table 1. 


was subsequently conducted in iqtree with default parameters. The four RDL 
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rates were calculated using the extant species richness of the associated family of 
each taxon, and the stem age (Ma) of that family taken from previously published 
time-calibrated phylogenies. Rdl copy number is provided within tip points 

and was determined by manual inspection of each genome and transcriptome 


Extended Data Fig. 2 | Species-level analysis of the relationship between Rdl 
duplications and net diversification rate. ASTRAL species tree of all insect 
species analyzed. The species tree was constructed using gene trees of 1123 
single-copy orthologs extracted using the BUSCO pipeline from published 
genome and transcriptome assemblies. Gene trees were inferred from aminoacid ^ assembly. 
alignments using maximum-likelihood estimation in IQ-Tree. Net diversification 
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Extended Data Fig. 3 | See next page for caption. 
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Extended Data Fig. 3 | View of the binding sites in the picrotoxin-bound 

RDL structures for each of the 2’ substitutions with the representative top 
poses (related to Fig. 1). (a to c) Model of the thymol (a-c)-bound Drosophila 
melanogaster RDL homo-pentamer. View of the binding pocket from the parallel 
tothe membrane plane (a). Side-on (b) and down-top (c) views of the picrotoxin 
and thymol-bound channel pore. (d) Plot of the pore radii in wild-type and 
mutant RDL models. A: wild-type. (e-1) The docking scores of the picrotoxin 


poses (top) are -7.09, -6.18, 77.73, and -7.18 in the structures of A2’S (e), A2’Q (f), 
A2'N (g), and AZ/P-T6'l (h), respectively. The docking scores of the thymol poses 
(top) are 5.02, —5.08, —5.25, and -5.12 in the structures of G279S (i), 1276F-G279S 
(j), V3391 (k), and A343T (D), respectively. Picrotoxinin, the main active ingredient 
of picrotoxin, and thymol are shown in ball-and-stick, amino acid side chains are 
shownassticks, and the dashed lines indicate hydrogen bonds. 
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Extended Data Fig. 4 | TraitRateProp analyses show the strong association between M22' residues of RDL and feeding states in insects. A rooted ultrametric 
phylogenetic tree, multiple alignments of RDL, and trait states from three orders (Hemiptera, Coleoptera, and Lepidoptera) were used, and see more information in 
method details. 
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Extended Data Fig. 5 | See next page for caption. 
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Extended Data Fig. 5 | Schematic of the strategy for generating knock-in 
alleles of Rdl using a two-step CRISPR-Cas9 gene editing method. (a) The 
genomic region of Rdl was replaced with a 1000-bp homology arm carrying a 
point mutation through homologous recombination. (b) The genomic region of 
Rdl was first replaced with a 3xP3-RFP marker. (c) The 3xP3-RFP marker was then 


replaced with ahomology arm carrying a point mutation. Yellow stars represent 
the point mutation. The crossing methods to generate knock-in lines are also 
shown. See also the generation of knock-in flies in Supplementary Tables 13-15 
for further details. Methods generally follow from refs. 41 and **. 


Nature Ecology & Evolution 


https://doi.org/10.1038/s41559-023-02127-4 


5 mM picrotoxin 


1» pato $5 ic 48 48 


Adult surivial (96) 


0 
0 12 24 36 48 0 


12 24 36 48 
Hours Hours 
[ed 5 mM thymol 10 mM thymol 


Adult surivial (96) 


0 T T T 1 0 T - ; 
0 12 24 36 48 0 12 24 36 48 
Hours Hours 
e 5 mM thymol 10 mM thymol 
S 
S 
E 
5 
o 
3 
"o 
< 
O 12 24 36 48 O 12 24 36 48 


Hours Hours 


10 mM picrotoxin 


te dedo tototo $m 


Io dodo dede do de 


ie to(81848 42 


Extended Data Fig. 6 | Amino acid replacements and duplications of Rdl 


show resistance to picrotoxin and thymol in gene-edited Drosophila 


melanogaster. (a-e) Adult survival of flies (see Methods for details) reared 


on diets containing picrotoxin (a and b) or thymol (c-e) among distinct lines. 
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The mutations in M2 (a-d) and M1 and M3 (e) were tested. 5-7 old virgin female 
flies were gently introduced into vials containing different concentrations 


of picrotoxin or thymol. Genotype: +/+, wild-type; S, A2/5; Q, A2/0; N, A2N;P, 


A2P; and I, T61. The log-rank (Mantel-Cox) test (Supplementary Table 16) was 
performed, mean + SEM (a, n = 6, 6,6, 6,6, 6, and 3 biological replicates; b, n - 3, 
3,3,3,3,3,3,and 3 biological replicates; c, n = 6, 6, 6,10,12, 5, and 5 biological 
replicates; d, n = 3, 3, 3, 3, 3, 3, 4, and 6 biological replicates; e, n= 3, 3,3, 5, 11, and 3 
biological replicates and 10 females per replicate). Experimental groups denoted 
by ‘ab’ are not significantly different from either a or b groups. 
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Extended Data Fig. 7 | The mutational path of Rdl in the lepidopteran lineages 
associated with increased net diversification of species increases resistance 
to picrotoxin in gene-edited Drosophila melanogaster. (a) Estimation of 
mean lifespan in adult female flies in CAFE assays at a range of picrotoxin 
concentrations. Means were compared using two-way ANOVA and post hoc 
Tukey's tests (Supplementary Table 16). (b) Estimation of LD50. Means were 


picrotoxin (mM) 


picrotoxin (mM) 

compared using two-way ANOVA and post hoc Tukey's tests (Supplementary 
Table 16). (c) Estimation of feeding rates. Changes were compared to the wild- 
typeflies. Means were compared using one-way ANOVA and post hoc Tukey's 
tests (Supplementary Table 16). Genotype: +/+, wild-type; S, A2/5; and Q, A2/Q. 
Each data point represents the mean + SEM of five biological replicates. *P « 0.05, 
**P < 0.01, ****P < 0.0001. 
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pentameric GABA receptors. Comparisons to wild-type flies are shown. n = 8,7, 
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Extended Data Fig. 9 | Daily activity across adult Drosophila melanogaster of a12-h light: dark cycle. Daily cross activity test in a 2-day window in the yellow 
different engineered genotypes. a-d, 5-7 old female flies of different genotypes box. mean + SEM, n = 24-32. Genotype: +/+, wild-type; S, A2/5; Q, A2'Q; N, A2'N; P, 
were gently introduced into each tube containing 2% agarose and 5% sucrose A2P;1, T61. 

diets, the other end was sealed with a cotton plug. Black and white bars represent 
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Extended Data Fig. 10 | The duplicated copies of Rdl are mapped on chromosomes of Acyrthosiphon pisum (Hemiptera) (a), Bombyx mori (Lepidoptera) (b), and 
Harmonia axyridis (Coleoptera) (c), respectively. Blue, green, and red circles represent Rdl1, Rdl2, and Rdl3, respectively. 
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